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Abstract 

The time-dependent density functional (TDDFT) equations may be written either for the Kohn- 
Sham orbitals in Hilbert space or for the single electron density matrix in Liouville space. A 
collective-oscillator, quasiparticle, representation of the density response of many-electron systems 
which explicitly reveals the relevant electronic coherence sizes is developed using the Liouville space 
representation of adiabatic TDDFT. Closed expressions for the nonlinear density-density response 
are derived, eliminating the need to solve nonlinear integral equations, as required in the Hilbert 
space formulation of the response. 
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I. INTRODUCTION 

Time dependent density functional theory offers a computationally tractable framework 
for calculating response functions of many-electron systems (such as molecules, semiconduc- 
tors and metals) to electrical or optical perturbations^ - — . The linear and nonlinear response 
functions (NRF) may then be used to extract useful information about electronic excited 
states (energies, transition and permanent dipole moments, etc.). The standard TDDFT 
algorithm involves the computation of the Kohn-Sham orbitals V'n(r) in Hilbert space,— and 
the response functions are calculated in two steps: First, a reference response is calculated 
for noninteracting electrons. Then a nonlinear integral equation is solved for the actual 
response function. The complexity of these calculations increases rapidly with the nonlinear 
order of the response. Alternatively, the Kohn-Sham equations may be written for the single 
electron density matrix p(r, r') and its evolution in Liouville space. We shall denote these 
as time dependent Hilbert space (TDHS) and time dependent Liouville space (TDLS) rep- 
resentations, respectively. The reduced single electron density matrix for N electron pairs 
ip n (r,t) is given by^ 

N 

p(r,r',t) = Y / Mr,t)r n (r',t) (1) 

?1=1 

A density matrix approach to the response was developed for the time dependent Hartree- 
Fock approximation (TDHF)^ - ^, which is formally very similar to TDDFT 2 -. A closed 
algebra of quasiparticles which satisfy a somewhat was established unusual scalar product, 
which allows a systematic order by order expansion of the density matrix in the fielcU^. In 
this article, we extend this approach to TDDFT and provide closed expressions for density- 
density response functions using the adiabatic exchange-correlation potential. 

Doubling the space dimensionality in Liouville space compared to Hilbert space (p(r, r') 
vs. ift n (r)) seems computationally prohibitive, in particular since within the Kohn-Sham 
approximation p(r,r') may be computed using ip n (r) using Eq.flTJ). However, recasting re- 
duced descriptions of many-body systems using density matrice o 15 ! 16 has many numerical 
advantages, as will be discussed later. In addition, the TDLS representation offers a valuable 
physical insight: (i) A quasiparticle picture of the response is naturally obtained in terms of 
the eigenstates of the linearized TDDFT. (ii) Closed expressions for response functions may 
be derived. In contrast, in the TDHS approach one needs to solve a chain of integral equa- 
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tions for each order self-consistently. (iii) The density matrix carriers explicit information 
about electronic coherence sizes. Molecular (e.g. Kohn-Sham) orbitals can be arbitrarily 
transformed to be localized or delocalized, without affecting the many-electron state. Their 
coherence length is therefore ill defined and does not carry a meaningful physical informa- 
tion. The density matrix is, on the other hand, uniquely defined. Optical excitations involve 
the creation of electron-hole pairs (excitons). The diagonal length of p(r, r') along r = r' 
reflects the center of mass motion of excitons whereas the off-diagonal length (along r — r') 
shows the size associated with the relative electronic-hole pairs motion. Plotting the density 
matrix provides physical insights and offers a real space N scaling description of electronic 
excitations when the density matrix is localized along r — r', which is typically the case™ 
iv) The response functions for any single-electron operator (e.g. current) can be readily 
obtained. This allows to compute, for example, conductivities. 

This paper is organized as follows: The time dependent Kohn-Sham equations of motion 
for the single electron density matrix are presented in Sec. II. The CEO (collective electron 
oscillator), quasiparticle, representation of the TDDFT density matrix is derived in Sec. III. 
The linear density-density response function is calculated in Sec. IV and the second-order and 
the third-order responses are given in Appendixes B and C, respectively. The advantages 
of the TDLS approach for computations of linear and non-linear time-dependent response 
functions, as well as possible extensions are summarized in Sec.V. 

II. TIME DEPENDENT KOHN-SHAM EQUATIONS FOR THE SINGLE ELEC- 
TRON DENSITY MATRIX 

The time dependent Kohn-Sham equations of motion for the density matrf^i are 
.dp(r,r',t) 



dt 

where we set (ft = 1) and 



L KS [p(r,r',t)](r,r',t) P (r,r',t); (2) 



p(r, t) = p(r, r, t) = n(r, t) (3) 

is the exact time-dependent charge density of electrons; L KS [p(r,r' ,t)](r,r f ,t) is the super- 
operator in Liouville space corresponding to the Kohn-Sham Hamiltonian. 

L KS [p(r, r', t)](r, r', t) = H KS [p(r, r, t)) (r, t) - H* KS [p{r', r', t)](r', t); (4) 
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Hks is the Kohn-Sham Hamiltonian 



H KS = H° KS +U 1 (r,t). (5) 

Hks i s the Hamiltonian of the isolated model and U\ represents the coupling to an external 
optical field. 

H° KS [p(v, r, t)](r, t) = -^ + U° KS [p(v, r, f)](r) + U (v); (6) 

m is an effective electron mass and the Kohn-Sham external potential is a functional of the 
charge density: 

U° KS [p(r,r,t)}(r) = J + f4 c [p(r, r, t)](r); (7) 

Z7o(r) is the field created by nuclei and U xc [p(r, r, t)](r) is the exchange-correlation potential 
in the adiabatic approximation. We further define the time-dependent external potential, 
given by U ext (r, t) = Uq(t) at time t <t and U ext (r, t) = Uo(r) + Ui(r, t) for t > to- 



III. QUASIPARTICLE REPRESENTATION OF DENSITY RESPONSE FUNC- 
TIONS 

The TDDFT equations describe the evolution of the reduced single-electron density ma- 
trix of a many-electron system driven by an external field. The calculation starts by comput- 
ing the ground state single electron density matrix p(r, r') in real space, which is determined 
by the stationary solution of the Kohn-Sham equation (Eq. (j2J). p is a key input to the 
TDDFT calculations. We then set p(r, r', t) = p(r, r') + Sp(r, r',t) where 5p(r,r',t) repre- 
sents the changes induced by the external field, described by the potential U\(r,t). The 
diagonal elements Sp(r, r, t) give the changes in charge density, whereas the off-diagonal el- 
ements p(r, r',t) represent the changes in electronic coherences between two points. Eq.© 
is then recast in the form 

i 96 ^^ = (^ 5 [p(r,r,t)](r)-^[p(r',r',t)](r') 

+ [/7 1 (r,t)-[/ 1 *(r',t)])p(r,r / ,t); (8) 

This equation may be solved for 5p(r, r', t) either in the frequency or in the time domain. 
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In order to calculate linear and nonlinear response functions we expand 5p(r, r', t) in 
powers of f/i(r,t): 8p(r,r f ,t) = Spi(r,r',t) + 5p2(r,r',t) + 5p 3 (r,r',£) + ... and solve the 
resulting equations of motion successively order by order. To that end we first expand Ug S 
in Taylor series in the vicinity of p to third order in Sp: 

U° KS [p(r,t)](r) = U° KS [p](r) +V f + V g + V h ; 



V f 



dr' 



r — r 



- + f xc [p]{v,v'))8p{v',t) ] 



V a = J dr'J dr"g xc [p](r,r',r")5p(r',t)8p(r",t); 

V h = J dr' J dr" J dr"'h xc [p](r,r',v",T"')5p(v',t)5p(r",t)5p(r"',t), (9) 



where for brevity we denote Sp(r, r, t) as 5p(r,t); e is the electron charge and f xc [p](r,r'), 
g xc [p](r, r', r") and h xc [p](r, r', r", r'") are the first, the second and the third order adiabatic 
exchange- correlation kernels^ (we consider the commonly used adiabatic approximation 
where they are assumed time-independent): 

su r . 



fxc[p}(r,r') 



5p(r',t') 



(10) 



r, r', r") 



S 2 U xc [p}(r) 



6p(r',t')5p(r",t") 



(11) 



/ap](r,r',r",r ,/A 



5 3 U xc [p](i 



5p(r',t')5p(r",t")Sp(r"',t"') 



We further define 



(12) 



f xc [p](v,v') = f xc [p](v,v') + 



r — r 



(13) 



and introduce the following notation for the product of two real space matrices C and r\ 



(C*J7)(r,r') = J C(r,r")r](r",r')dr". 
For an operator A[£(r, r)](r) and a matrix ^(r, r') we denote 

[A, r)](r, r') = (A[C(r, r)](r) - A*[C(r', r')](r')) v (t, r'), 
and the commutator of two single-electron matrixes ( and rj 



(14) 



(15) 



[C ) r?](r,rO = (C*?7)(r,r / )-(?7*C)(r,r / ). 



(16) 



Since p(r,r',t) represents a single Slater determinate, it must be idempotent, and, conse- 
quently, not all elements of 5p(r,r',t) are independent^. When expanded using the ground 
state Kohn-Sham orbitals, Sp can be divided into four blocks: electron-hole and hole-electron 
(interband) and electron-electron and hole- hole (intraband). Only the interband elements 
of £(r, r', t) are independent and the intraband blocks T can be expressed in terms of these 
element o 24 i 25 i 26 . We thus decompose Sp(r,r',t) as (see Appendix A) 

6p(r,r',t) = ary,t) + T{£{r,r',t)), (17) 

It will be desirable to remain in real space and avoid switching back and forth to the orbital 
basis. Fortunately, this goal can be accomplished and £ can be computed directly in real 
space using the double commutator: 

ttr,r',t) = [[5p(t),p),p)(T,r'). (18) 

The operation in the r.h.s. is a projection operator that projects Sp into the interband 
subspace (see Eq. fjAlj) ) 18 . Once £ is calculated we obtain for T (see Appendix A) 

T(£(r, r ', *)) = (/ - 2p) * (£(t) * £(*) + £(t) * £(t) * * £(t) + ■ ■ -)(r, r'), (19) 

where all £(r, r', t) are taken at time t. 

The next step is to develop a convenient quasiparticle algebra for the independent dy- 
namical variables £ which represent collective electronic excitations. To that end we start 
with the eigenvalue equation corresponding to the linearized Kohn-Sham equations 

LUr,r') = n a Ur,r), (20) 

where 

L = H f (r)-H* f (v'), (21) 

H f (r)Ury) = - V ^ ,r0 + (/ dv'J dv" f xc [^ OUr', r")) ?(r,0, (22) 

and 

/L[p](r,r' ; r") = /L[p](r,r')5(r'-r"). (23) 
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The eigenmodes £ Q come in pairs. Each pair of conjugated modes £ a (r, r') and ^(r, r') 
forms a collective electronic oscillator a (where a = 1, 2, . . .). We define £- Q (r, r') = ££(r, r') 
and = — Q a where a = ±1, ±2, . . .. The algebraic properties of these oscillators were 
developed in Ref.— and reviewed in Ref.— . 

We further introduce the following scalar product of any two interband matrices £ and r\. 

(fa) = J ' dvj dv'p[^,r])(r,r')S(r-r'). (24) 

The bra (ket) notation underscores the similarity with Dirac's Hilbert space notation. 
Eq. (J21)) satisfies the following relations: 

(£fo> = { V W = -(v\0 . (25) 

Using this scalar product, we shall normalize these modes as follows 

(&\Zf>)=6ap, (U^) = 0. (26) 

Eqs. ()25|) and (|26|) show that this is an unusual scalar product: for one thing, the norm of a 
vector is zero. Nevertheless, it has one important property: it allows us to expand £(r, r', t) 
in terms of collective electron oscillator (CEO) modes £ Q (r, r') 

£(r,r',*)= ]T £ a (r,r'> a (t), (27) 

a=±l,±2... 

where z a (t) = an d its complex conjugate z.. a (t) = z^(t) constitute complex os- 

cillator amplitudes. £ Q thus form a complete basis for representing an arbitrary interband 
matrix. 

Combining Eq.()27|) and Eq. (fT7|) . we get 

6p(r,r',t) = ^p a (Y,r')z a (t) + ]-^p a p(r,Y')z a (t)z^(t) 
+ 3 ^2pa,p^,r')z a {t)z,3{t)z^{t), 

a/37 

a,/?, 7 = ±l,±2,... (28) 

where we only retained terms that contribute to the third order response. The expansion 
coefficients in the r.h.s. are given by 

p a (r,r') = £ a (r,r'), (29) 
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p a p(r, r') = (/ - 2p) * (£ Q * £p + £p * £«)(r, r'), 



(30) 



Pa,/3 7 (r, r) = -f a * (£g * £ 7 + £ 7 * £s)(r, r'). (31) 

The amplitudes of the adjoint (negative frequency) variables are simply the complex conju- 
gates (z„ a = z* a ). 

Upon the substitution of Eqs. (|2B |) .(|25 j) .(j3D )l . and (j3*Tj) into Eq. (|SJ) we obtain the equation 
of motion for the CEO amplitude z(t) to third order in the external field: 

.dz a (t) 



ot 



where 



Q a z a (t) + K- a + K_ a p zp{i) + Y K^ af3l zp(t)z 7 (t) + ^ ^-q/? 7 5 zp{t)z~,{t)zs(t)\ 

P P 1 p-y5 

a,/3,7,5 = ±l,±2,..., (32) 



K- a = J U 1 (r,t)p^ a (r)dr; P- a (r) = p_ Q (r,r); 

K- a p = J Ui(r, t)p_ Qj/3 (r)dr; p- a ,p( r ) = P-a,p( r , r ); 

K_ a/3y = Uicsi-cPi) + / ^i( r > t)p- aj p 7 (r)dr; p~ a ,p y (r) = p- a> p-y(r, r); 

Uk S is the Kohn-Sham external potential, expanded to the third-order kernel Eq. Q, and 
U KS(- a py) anci U KS(- a p-y6) are given in Eq. (jB6j) and Eq. (JOSJ), respectively. 

Eqs. (|3~2|l together with (j2Hj) constitute the quasiparticle algorithm for computing density 
response functions. These nonlinear equations which map the system onto coupled collective 
classical oscillators, may be solved by expanding z(t) (z*(t)) in powers of the external field 
E/iM): 

z (t) = (t) + z i2) (t) + z® (t) + ..., (34) 

Substituting this in Eq. ()32j) we can successively solve for z^ order by order which when 
substituted in Eq (|28J) will give Spj. The present truncation of Eqs. (}28|) and (|32J) allows to 
compute response functions up to third order. However, this approach may be extended to 
NRF of arbitrary order by simply adding higher terms to Eqs. (|2*K|) and (|32|1 . 
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IV. THE LINEAR DENSITY RESPONSE 

The following calculation of the linear response illustrates the strategy for computing 
response functions. The first-order density matrix (Eq. (jSJ)) satisfies: 

t d5pl{r dt r '^ = L5 Pl (r,v',t) + [^(r.f) - U* 1 (r',t)]5p 1 (r,r',t). (35) 

The equation of motion for z^>{t) is obtained from Eq. (jHSJ) together with Eq. lfTTj) and the 
expansion Eq. (fTT)|) : 

i^Of 1 = V a z i a ) (t) + j U 1 (r,t)p- a (r)dr, a = ±l,±2,... (36) 

The solution of this equation can be represented as 

4 1} (*) = is « f dT I drU 1 (r,r)p_ a (r)G a (t-T), (37) 

J — oo J 

where positive and negative a correspond to z^\t) and z^\t), respectively, s a = sign(a), 
and we have introduced the time-domain Green function 

G a (t) = e{t)e~ inat , G- a {t) = 9{t)e- in - at = 6{t)e inat , (38) 

where 9(t) is the Heavyside step function. 

Inserting Eq. (JBTj) into Eq. (|2~%j) we obtain for the density to linear order 

6 Pl (T,t) = 5 Pl (T,T,t)= Yl z( a\t)Pa(r). (39) 

a=±l,±2,... 

This can be recast in the form 

S Pl (r,t) = f dr [dT , U 1 (T>,T) X ( 1) (t,T,Ty), (40) 

J — oo J 

where the first order time-domain linear density-density response function is given byii: 

X (1) {t,T,r,r') = i is aP -a{r') P a{r)G a (t-T). (41) 

a=±l,±2,... 

The corresponding frequency- domain density- density linear response function 
X^(— ^s', is defined by 

/oo fj/.j r 
— dv'x^H-^u^v'^ir'^). (42) 
-oo /7T J 



Here U\(r,u) is the Fourier transform of the time-dependent external field Ui(r,t) defined 

as 

/(r, u) = J dtf(r, t)e™ , /(r, t) = ^ J dwf{v, u)e~ wt . (43) 
By comparing Eqs. (j^Ojl with Eqs. (|42jl and using Eq. we have: 

rite**' / dre-^x {1 \t,T,r,r'), (44) 

-oo J — oo 

The linear order response functions is usually denoted^ 

X (1) K = w;w,r,rO = 27r5(-^ + a;)x (1) (^,r,r'), (45) 
Using Eqs. (}1H) and we obtain the linear density-density response function 

x a Wy) = E y-.(^(r) = E 2^(W _ (46) 

a=±l,±2,... a W af=l,2,... W 

Here and below Q a is positive (negative) for all a > (a < 0), following the convention 

Finally, the static linear density- density response can be obtained from Eq. (j4T)j) by setting 
uo = 0: 

X (1) (0,r,r') = £ W^KCO (47) 

a=l,2,... 

Higher order response functions may be computed in the same way and depend on 
the same ingredients: the collective electronic oscillator (CEO) modes. Closed expres- 
sions for the second-order density-density response functions are given in Appendix B 
(Eqs. (|B12|) . (jB19|) and (|B20|) ) and the third order response is presented in Appendix C 
(Eqs. (IHTol) . (IC231) and dC32|) V 



V. DISCUSSION 



The density matrix TDLS representation of linear and non-linear TDDFT density-density 
response functions has several notable computational advantages compared with the Hilbert 
space (orbital) TDHS representation. To facilitate the comparison, we have outlined the 
TDHS in Appendix D. In the TDLS approach we need to solve the CEO eigenfunction 
and eigenvalue problem in Liouville space (Eq. (j20J0 to obtain the quasiparticle spectrum 
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of excitations (the CEO modes)^. Eqs.©, flEEl and 

can then be directly used 

to calculate time-dependent response functions in a single step. In TDHS, in contrast, we 
proceed in two steps: we first solve K Kohn-Sham equations Eq. (jDl|) with Eq. (|D2|) for the 
Kohn-Sham orbitals, and calculate the linear response for non-interacting particles Eq. (jD3|) 
and as a second step, we solve the Dyson-type integral equation Eq. (jD4|) to obtain linear 
response 11 . Calculating the non-linear response in this approach is even more complicated. 
For example, obtaining the second-order response requires solving the Dyson-type integral 
equation Eq. (|D7[) using the integral equation for the linear response Eq.(jD6]). Similarly, 
obtaining the third-order response requires the solution of Dyson-type integral equation 
Eq. (|Dll|) using the integral equations for the second Eq. (jD10|) and linear response Eq. (jD6|) . 
The computation of non-linear higher-order response thus involves the chain of self-consistent 
integral equations for all previous responses, while the TDLS provides a closed expression 
for the response functions. Clearly, the TDLS approach provides a much faster one-step 
algorithm for computations of response functions. 

The CEO modes are nonlocal and act in a six dimensional space (r, r') compared to the 
three dimensional space (r) of the Kohn-Sham orbitals. However, several points help reduce 
the computational cost. First, we can reduce the number of Liouville space equations by 
solving only for particle-hole (interband) density matrix £ instead of the entire density matrix 
5p: the particle-particle and hole-hole (intraband) density matrix T(£) can be obtained from 
Eq. (jl9j) . which follows from the idempotent property of single Slater determinant. For a 
basis set of K orbitals the number of elements is the KN rather than K 2 . Moreover, the 
density matrices have non- vanishing elements only when |r — r'| is less than a coherence size, 
which is typically very short. This allows to neglect many density matrix elements, further 
reducing the size. It is not possible to include the coherent size in TDHS computations. 

The present work can be extended in several ways. Even though we have derived closed 
expressions for the response up to third order, the computation of higher-order responses is 
straightforward and merely involves carrying out the expansions in Eqs. (J2%|) and (|H2j) . The 
procedure, used here to compute response functions of the charge density, can be used to 
obtain response functions for any single-electron operator. For example, conductivities are 
given by the current response 

f\r,t)=Tr((v6p^(t))(r,r')), (48) 
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where v(r,r') = — 2^;(V r ' — V r ) is the velocity operator. 

Another interesting possible extension of this work is to include non-adiabatic exchange- 
correlation potentials, as outlined recently for the linear response 2 ^. In general, the exchange- 
correlation potential and exchange-correlation kernels are time-dependent 2 ^*. This time de- 
pendence has been neglected within the adiabatic approximation used here. If we relax this 
approximation, the eigenvalue equation for the Liouville superoperator L, Eq. (J20J), should 
be replaced by2i 

L(Q a )Ur,r') = Q a Ur,r'). (49) 

Methods for solving Eq. (j49j) using the frequency-dependent functional of Gross and Kohn^ 
were described in^l. 

The time dependent Hartree-Fock TDHF (also known as RPA) is an alternative widely 
used approximation for computing response function a 24 i 25 i 26 . Within the density matrix 
representation the TDHF is formally equivalent to TDDFT, provided we use the following 
exchange-correlation potential that depends on the density matrix (rather than merely on 
the charge density)^: 

1 5 

Vks-tdhfMt) = - 2 j^ 

By substituting Ux S [p](r) = Ux S _ TDHF [p](r) from Eq. (|5Dj) into Eq. (JJJJ), we obtain the 
TDHF density-density responses. The TDLS formalism and the associated CEO algebra 
were first derived for the TDH F 18 i 20 and this paper extends these results to TDDFT. The 
difference between TDDFT and TDHF equations is that in the former each order of the 
perturbed density matrix has its own equation of motion with a different electron-electron 
interaction potential, because the TDDFT exchange-correlation potential, expanded in terms 
of exchange-correlation kernels, is different in each order, while in TDHF the electron- 
electron interaction is the same for all orders 20 . 
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dr / dr 



/P(r,r>(r',r) 



(50) 



n(r') = o(r / .r f ) 
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APPENDIX A: PROJECTION OF THE DENSITY MATRIX INTO THE INTER- 
BAND SUBSPACE 

Since the many-electron wave function is represented at all times by a single Slater de- 
terminant, the total density matrix p(r, r', t) must be a projector. The idempotent property 
((p) 2 (r,r') = p(r,r')) allows us to project any single particle matrix ( into the interband 
(p-h) subspace^ 

C P -,(r,r') = [[C,p],p](r,r'). (Al) 

Consequently, not all elements of the density matrix are independent^ 8 -. The idempotent 
property gives 

((p + Sp(t)) * (p + 5p(t)))(r, r') = p(r, r') + 5p(r, r', t). (A2) 

The number of degrees of freedom of 5p subject to the condition Eq. (|A2|) is precisely the 
number of its particle-hole matrix elements, and T(£(r, r', £)) can therefore be expressed in 
terms of £(r, r',t). Using Eq. (jA2|) and Eq. (fTTj) we get 

(p + £(f) + T(t{t))) * (p + m + T(m))(r, = P(r, r') + £(r, r', t) + T{t(r, r', t)). (A3) 

To simplify this expression we use the following relations, which follow from Eq. (jAl|) : 
p(r, r') = (p*p)(r, r'), £(r, r', t) = p*£(t)(r, r')+£(t)*p(r, r') and T(t)*p(r, r') = p*T(t)(r, r'). 
The following rule may be applied to separate the remaining terms: product of two inter- (or 
two intra-) band matrices gives an intraband matrix, whereas product of inter- into intra- 
(or intra- into inter-) band matrices results in an interband matrix. Finally, the intraband 
part of Eq. (JX3J) is 

(Tim) * T(£(f)))(r, r') + (2p - /) * T(£(*))(r, r') + (£(*) * £(t))(r, r') = 0, (A4) 

where I is the unit matrix in the real space (I(r, r') = 5(r — r')). The formal solution of this 
quadratic equation, with the condition T(£(r, r', £) = 0) = yields 

T(£(r,r',t)) = (p - * (/ - ^ - 4£(t) * ^)) (r,r'). (A5) 

Eq. (fT^j) is obtained by expanding Eq. (|A5|) in powers of £. 
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APPENDIX B: THE SECOND-ORDER DENSITY RESPONSE 



Expanding Eq. (JHJ) to second order, we obtain 
.<9<5p 2 (r,r',t) 



i- 



L5p 2 (r,r',t) + (V f [S Pl ](r,t) - V f *[5 Pl ](r', t))6 Pl (v, r', t) 



+ (V;[S Pl S Pl ](r,t) - Vg*[5 Pl 6 Pl )(r',t)) P (r,r') + {U^t) - U{(r',t))5 Pl {r,r',t), (Bl) 

where 

V f [5 Pl ](r,t) = J dr"f! xc [ P }(r,r , ,r")5p l (r',r",t), (B2) 

and 

V$ Pl 6 Pl )(r,t) = J dr'J dv" j dv>" J dr""g xc [p](r, r', r", r'", r"") 

5 Pl (r',r",t)5 Pl (r w ,r"",t), (B3) 

where 

g xc [p)(r, r', r", r'", v"") = g xc [p](r, r', r>")8(r> - r")5(v"' - v""). (B4) 
Using Eq. (J32J) . the equation of motion for is 



n«4 2) (*)+£4 1} (*) / ^(r,t)p-^(r)rfr 



at , 

/3 7 

a,/3, 7 = ±l,±2,..., (B5) 

where 

K(a^y) = ^Tr((/-2p)*((e /3 *e 7 + e 7 *^)^(U 

+ (&> * & + & * £a)^(£ 7 ) + (a * £ 7 + £ 7 * Z a )V g (£p))\ (B6) 

where 

= V,fcJ(r)& + VJIUJWP, (B7) 

and 

Tr(C(r,r')) = J dr J dr'((r,r')5(r - r') = / C(r,r)dr = / COO*"- ( B 8) 
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The solution of Eq. (jB5|) is 

4 2) (*) = * f dr x JdrJ dr'saG^t-n^in^r'), (B9) 

where 

r^(ri,r,r') = J2 V 9(~<*Pi) I dr i I d nUi{r, ri)C/i(r', r 2 )p_ /3 (r)p_ 7 (r / )s /3 s 7 
Gpin - r 2 )G 7 (ri - r 3 ) + ^ it/i(r, ri)p_ a/3 (r) 

f 1 dr 2 U 1 (r',T 2 )p-f ) (r')8f i Gf i {Ti-T 2 ). (BIO) 

J— oo 

The time domain second-order density-density response function x {t> r i; r 2) is defined^: 
8p 2 (r,t) = 5p 2 (r,r,t) = \ f d n f dr 2 f dr' f dr"?7 1 (r , ,r 1 )?7 1 (r",r 2 ) 

Z J —oo J —oo J J 

X {2) (t,T U T 2 ,r,r',r"). (Bll) 

Inserting Eqs. (J37)l and (jB9|) into Eq. (J28j) and keeping all terms up to the second order we 
find that the second-order response function finally has three contributions: 

X {2) (t, n, r 2 , r, r', r") = X f\t, r lt r 2 , r, r', r") + n , r 2) r, r', r") 

+ X?}i(-t,Ti,T 2 ,T,v',r"), (B12) 

where 

^(t.n.Ti.ry,^) = 2^p_ Qi/3 (r)p Q (r')p- /3 (r") 

a/3 

s a spG a {t - ti)G p {tx - r 2 ), (B13) 
dr V 9(-^)Pa(r)p^f3{r')p^(r") 

s a Sps 7 G a (t - r)Gp(r - r 1 )G 7 (r - r 2 ), (B14) 
xSr^ri.^ry.r") = 2 $> Q/3 (r)p_ a (r')p-Mr") 

a/3 

s a Sf,G a (t-Ti)Gfi{t-T2). (B15) 

The corresponding frequency domain density- density response function x^ 2 \~ ^s] k>i, ^2) 
is defined by 

1 Z" 00 GL>i Z" 00 



fo(r, W .) = \f_J^ l~J^ J dr'J dv" X ^\-u s ,u l ^vyy) 

fMr'^OCMr",^). (B16) 
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The relation between response functions and charge densities are obtained by comparing 
Eq. ()B11|) with Eq. (jB16|) and using the Fourier transform Eq. (jmj) : 



X (2) (-^;^i,^2,r,r',r") 



(°° dte iUat f dr ie -^ lTl f dT 2 e- iu,2T2 x {2) (t,T 1 ,T2,r,r , ,r"). (B17) 

J —oo J— oo J —oo 

The second order response function is usually denoted 

X^i^s = ui + u 2 ,Ui,U2,r,r',r") = 2ttS(-u s + oj\ + u) 2 ) 

X^K^ry,!*). (B18) 

Using Eq. ()B17|) we finally obtain the second-order density- density response function 
which is symmetric with respect to U\ and u 2 permutations 

V (2)A, r / _ o V ^(- a /37)P«( r )P-/3( r >-7( r "K^ 

X v ^i,^ 2 ,r,r , r ) - A 2^ 7?; \7?S 77n V 

+ P- a/3 (r)p Q (r')p_ /3 (r")s a s / 3 

a/3 

+ E 

a/3 





- ^2)(fy3 


- wi) 


P-a/3(r)Po 


,(r')P-/3(r' 




- Ui 




- W2) 


Pa/3 (r) Pa ( 




SaS/3 



Here and below fij,, i/ = a, f3 is positive (negative) for all v > {y < 0) following to the 
convention Q- u = — Q u . 

By setting uj\ and uj 2 to zero and using the identities s u fl u = 1^1 and p_ u = p u we obtain 
the second order static density-density response: 

X (2) (0 r r' r") = 2 V ^W^P^WP/^'Wr") , 3 ^ P^Wpo(OMO 

a,/3,7 = ±1,±2,.... (B20) 



APPENDIX C: THE THIRD-ORDER DENSITY RESPONSE 

Expanding Eq. (jHJ) to third order, we obtain 

+ (^p 2 ](r,t) - V7[5p 2 ](r',t))5 Pl (r,r',t) 
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+ (V^[8p 1 8p 1 )(r,t) - V^[5p 1 5p 1 ](r',t))6 Pl (r,r',t) 

+ (V;[5 Pl , 8p 2 )(v, t) - V' g *[8 Pl , 8p 2 ](r', t) + V^[8p 2 , 8 Pl )(r, t) 

- Vg*[5p 2 ,5 Pl }(r',t))p(r,r') + (V h [6p 1 ,6p 1 ,Sp 1 }(r,t) 

- V:[5p 1 ,5p 1 ,8p 1 ](r',t))p(r,r') + (^(r,t) - U^(r',t))8p 3 (r,r',t), (CI) 

where 

V h [8p 1 8p 1 8p 1 )(v,t) = J dr'J dv" J dv'" j dv"" J dv'"" J dv"""h xc [p)(r, r', r", r'", v"", r"'", r""") 

5 Pl (r', r", t)<5 Pl (r"\ r"", ^ Pl (r'"", r""", f), (C2) 

where 

^ c [p](r, r', r", r'", r"", r'"", r""") = h xc [p}(*, r', r'", v'"")5(v' - v")5(v'" - v"")5(v""' - r"""\C3) 
Using Eq. the equation of motion for is 



dt 



= Sl a z®(t)+zf\t)Y, / U 1 (v,t)p_ a ^v)dv 
+ E^W^W / C/i(r,t)p_ a/37 (r)dr 

(37 

+ 2^^ ( _^ 7) 4 1) (t)4 2 )(t) 

/3 7 5 

a,/?,7,5 = ±l,±2,..., (C4) 



where 



V^y*) = -Tr((7 - 2p) * * 6 + 6 * - 2p) * * £ 7 (C5) 

+ ?7 * &))) + " 2p) * (£„ + 6») W - 2p) * (6 * ^ 

+ ^ * &))) + ^Tr((J - 2p) * (£ Q * & + C/9 * &*) W - 2p) * (& * &* + £ 7 * 6))) 
o 

- l -Tr({t a V h (^) + * (£ 7 * is + 6 * £y)) 

- ^Tr((e Q v,(e 7 ) + * * & + 6 * &)) 

- ^Tr((£a^(&) + * * £ 7 + £ 7 * &)), (C6) 
The solution of Eq. (|C4|) is 

(t) = i /' rf Tl / drsAit-r^i^r), (C7) 

J —00 J 
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where 





Pi 

+ EVn-a^ir^in^in), (C8) 

/3 7 <5 

where z^(ti) and z' 2 )(ti) are given by Eqs. (j37p and ()B9|) . 

Time domain third-order density-density response function x^(£, ri, . . . , T3) is defined^: 

Sp 3 (r,t) = Sp 3 (r,r,t) = ~ f d n f dr 2 f dr 3 f dr' f dr" f dv'"U x {v' , n) 

O J— 00 J -co J— 00 J J J 

fMr", r 2 )[/ 1 (r'", r 3 ) X (3) (^ n, 75, r 3 , r, r', r", r'"). (G9) 

Inserting Eqs. (}3T|) and (jB9|) and ()C7|) into Eq. (J2Hj) and retaining all terms to third order 
we obtain the final expression for the third-order response function: 



X ^(t, Tl , r 2 , r 3 , r, r', r", r'") = X f + X f} + X f n + X fj + X f + x §\ + X % + X $ n , (CIO) 
where 

xi 3) (^ r i, t 2 ,t 3 , r,r',r",r"') = -6z ^ p_ a/3 (r)p_ /37 (r / )p a (r") 

a/37 

p_ 7 (r'")s a s /3 s 7 G' a (t - Tx)Gp{Ti - t 2 )G 1 (t 2 - t 3 ), (Cll) 

X//(^ r ir2r3,r,r',r",r w ) = 12 ^ P-a /3 (r)V 3( _ / 3 7(5 )p Q (r')p_ 7 (r")p-5(r" / )s Q s / 3S 7 s ( 5 

a/375 

x f t drG a (t-T 1 )Gp(r 1 -r)G 1 (T-T 2 )G s (T-T 3 ), (C12) 

xS7(^ r i,T-2,7-3,r,r',r",r w ) = -62 ^ p_ Q07 (r)p a (r , )p-^(r") 

a/37 

p^{r"')s a SpS^G a {t - ri)G /3 (ri - r 2 )G 7 (r 1 - r 3 ), (C13) 

X/y(^n,r2,r3,r,r',r",r w ) = 12^ V 3( _ a/37 )p_ 75 (r)p a (r / )p- / 3(r / ')p-5(r" / )s a s /3 s 7 s 5 

a/375 

x /' < drG a (t-r)G /3 (r-r 1 )G' 7 (T-r 2 )G' 5 (r 2 -T3), (C14) 
Xy ) (t,r 1 ,r 2 ,T 3 ,r,r',r",r'") = 12i ^ K, ( _ a/ 3 7) K, ( _ 7 ^ ) p a (r)p_ /3 (r')p-5(r") 

a/375r) 
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x f dr f T dr'G a {t - T)Gp{r - r 1 )G 7 (r - r') 

G 5 (r'-r 2 )G ri (r'-r 3 ), (C16) 

Xv ) I (t,T U T 2 ,T 3 , r,r',r",r"') = 6 V^ ( _ Q/375) p Q (r)p_ /3 (r')p^ 7 (r")p_ 5 (r w )s Q s /3 s 7 s ( 5 

a/375 

x f dTG a (t-T)Gp(T-n)G 7 (T-T 2 )G s (T-T 3 ), (C17) 
Xy!/(t,ri,r 2 ,T 3 ,r,r',r",r w ) = -12i ^ p a p{v)p_^{v')p_ a {v") 

p_ 7 (r" / )s a s /3 s 7 G' a (t - Ti)Gp{r - t 2 )G 7 (t 2 - r 3 ), (C18) 



xS/ 3 m(^i,r2,T 3 ,r,r',r",r'") = 12 51 Pa/3(r)V 3 (_ /37 5)p_ a (r')p 7 (r")p-5(r / ")saS/3S 7 s ( 5 



x 



a/375 



x / drG a (t-T 1 )G (t-r)G 7 (r-r 2 )G s (T-r 3 ). (C19) 

The corresponding frequency domain density-density response function 
X^ 3) {-uj s ]uji, ...,u 3 ) is defined by 

1 f 00 duoi f 00 duo 2 f°° duo 3 



O Z7T J-00 2% J-00 Ltx J J J 

U 1 (r',u 1 )U 1 (r",co 2 )U 1 (r'",co 3 ). (C20) 

The relation between response functions and charge densities are obtained by comparing 
Eq. (jC9|) with Eqs. (jC20j) and using the Fourier transform Eq. (^3]) : 



X {3) {-u s ; u)i,u) 2 ,u) 3 ,v,r',r", r'") = 

/oo 
dfe w / dr ie - luJ1Tl I dr 2 e' luJ2T2 I dr 3 e 
-OO J — OO J -co J —00 

X W{t,n,T 2 ,T 3 ,r,r'yy). (C21) 



OO r-t rt r-t 

iu B t I twin / r /, r „ a - iu >2T2 / ^„ C -*^3T'3 



The third order response function is usually denoted 

X (3) (^ s = cJi +o; 2 + u) 3 ,uJ u u) 2 ,uJ 3 ,r,r\r",r'") = 2tt5(-uj s + u x + lu 2 + u 3 ) 

x^i^^^ryyy). (022) 

Using Eqs. (|C21|) and (|C22|) we obtain the following 8-term expression for the third-order 
density response function (symmetrized with respect to Ui,u 2 , and u 3 permutations) 

perm 

X (3) (^i,^^3,r,r',r",r'")= £ (xf + xff + X?A + ■ ■ • Xvm) , (C23) 

19 



where 

(3) = y P- a p(Y)p- l 3 1 (Y')p a {r")p- 1 {v'")s a SpS 1 ^ C24 ^ 

^(3) = V- P- Q /3 (r) K,(-/3 7 ^) Pa (r>-7 (Qp-j (r W ) S a SpS y S S ^ 

(3) = P- Q/?7 (r)p a (r>_^(r // )p„ 7 (r w )g a g /3 g 7 ^ C26 ^ 

Y3) _ 2yg ( „^ 7 )p_ 7 ^r)p a (r / )p_ /3 (r // )p_^(r w ) SgSps^ss (C27) 

a/3 7 <5 (^« - ^1 - ^2 - ^3) (^/3 - Wl) (^7 - a> 2 - ^3) (^<5 ~ ^3) ' 

(3) 

Xv = 

y 2y g( _ a/37) y g( _ 7 ^ ) p a (r)p_ /3 (r / )p_g(r // )p_ r? (r w )g a 5 /3 g 7 g^ q 

Jh6 V - ^1 - ^2 - ^3)(fy3 - ^l)(^ 7 - ^2 - ^3)(^5 - ^2)(^r, - W 3 ) 

(3) = Vft.(- a ^)P a (r)p- / 3(r / )p_ 7 (r // )p_Xr /// )g a g /3 g 7 g^ 2g 

^ c^(J (^« - ^1 - ^2 - ^3) (ty(3 - Wi) (f) 7 - U 2 ) (Q S - U 3 ) ' ~ 

(3) = v Pa/3 (r) p_/3 7 (rQp-g (r") p_ 7 (r'") SgS^ 
Xvn (O a - a*) (fy - a* - a*) (O, - ws) ' 1 ) 

(3) = y Pa/3 (r) ^(-/frtf) P-a (Q p_ 7 (r") p_j (r w ) SgSpS^ss 31 

Here v = a,/?, 7,5,77 = ±1,±2, ... and Vt v is positive (negative) for all v > (y < 0) 
according to the convention Q_ v = —Q u . Note, that in Eq. (jB19|) the permutations over u>\ 
and u 2 were written explicitly. Finally, by setting u%, uj 2 and o> 3 to zero and using identities 
s v VL v = \VL v \ and p- v = p u we obtain the third order static density- density response: 

X {3) (0 r r' r" r'") = 6^ P^Mp^')p^")Pi^'") , g y 2p^(r)p^ 7 (r')p a (r")p 7 (r w ) 



a/3-y l"""P°°7l a/3 7 



+ 6 y 4p Q/ 3(r)V A g( _^)p a (r / )p 7 (r // )pX 
+ g y» 2F g(a/37) V^-^) p a {v)pp{T')p s (r")p v (r"') 

+ g Vft(a^) gg ( r ) Pj ( r ^7 ( r// ) gg ( rW ) 

a,/5, 7 ,5,r/ = ±l,±2,.... (C32) 
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APPENDIX D: THE HILBERT SPACE (TDHS) REPRESENTATION OF RE- 
SPONSE FUNCTIONS 



The ordinary ground-state Kohn-Sham equations in Hilbert space ar o 2 i n (h = 1) 

#KsM r )](r)y?;(r) = ^-(r); (Dl) 

D%W| 2 = Mr); (D2) 
j=i 

where n (r) is a true initial ground-state charge density of electrons (the charge of an electron 
e = 1; HKs[ n o( r )]( r ) is the Kohn-Sham Hamiltonian Eq.fJUJ) with the time-independent 
external field U ext (r) = Uo(r). 

The density-density linear response Eqs. (140)1 . (jBllj) . ()C9)) in the TDHS framework can be 
calculated in two steps, using a perturbative expansion in the electron-electron interaction^. 
The first step is the calculation of the response function x^(a;,r,r') of non-interacting 
particles with unperturbed density n Q in terms of the static unperturbed Kohn-Sham orbitals 

M r ) 

X U(u,r,r>) = , _ t _ ■ (D3) 

j,k \ L j tk > * t 'l 

Here, {fk,fj) are the occupation numbers (0 or 1) of the KS orbitals; Sj and are corre- 
sponding energy levels of the non- interacting particles. The summation in Eq. (JD3|) extends 
over both occupied and unoccupied orbitals, including the continuum states. 

The linear response function for the system of interacting particles x {u, r, r') is obtained 
in a second step by solving Dyson-type integral equation: 



X (1) (^,r,r') = x. 



i 1 )(^,r,r') + jdxjdx' X «(^r,x) 

7| + /kN(w,x 1 x / )| X (1) (^,x',r'), (D4) 



X — X' 



where f xc [no](u, r,r') is a Fourier transform with respect to time of the time- dependent 
exchange-correlation kernel f xc [n }(r,t,r' ,t'): 

5U xc [n](r,t) 



f xc [n ](r,t,r ,t 



5n(r',t') 



no 



(D5) 



21 



and the exact frequency-dependent linear density response is given by 



ni(u,r) = J dyxi 1) (w,r,y)f/i(y,a;) 

+ J dy J dy' X ^(u,r,y) (r^-^ + fxc[n }(u, y, y'^jn^u, y')- (D6) 

The Bethe-Salpeter-type integral equation for the second order density- density response is 
represented in four- dimensional coordinate-time space 



X {2) (x,y,y') = JdzJ dz' X f\x,z,z>) 



SU KS (z') 



5U ext (y>) 

+ j dzx^(x,z) j dz'j dz"g xc (z,z',z")x^(z',y) X ^(z",y') 

+ ( dz X ^(x,z) [ dz'(w(z,z f ) + fU*,*))x®W,V,l/), (D7) 



;ds) 



where the time-dependent second-order exchange-correlation kernel g xc is defined as: 

5 2 U„Jz) 



g xc (z, z, z") 



no 



Sn(z')Sn(z") 

the kernel of the (instantaneous) Coulomb interaction w(x, x') is 

e 2 5(t - 1') 



w(x, x') 



r — r 



(D9) 



The time-dependent Kohn-Sham equations for the second-order density response are 
finally obtained by inserting Eq. (|D7|) into Eq. (jBll|) : 

n 2 {x) = - J dz J dz'xf > {x,z,z , )U K s[ni]{z)U KS [n 1 ]{z') 

+ ~ J dz J dz' J dz" X { }\x,z)g xc {z,z\z")ni(z')ni{z") 

+ JdzJ dz' X V( Xl z)(w(z,z') + f xc {z,z'))n 2 {z'). (D10) 

Solving Eq. (|D6j) first, allows for the subsequent solution of the self consistent Eq. (|D10J ) 
which is quadratic in the (effective) perturbing potential UKs[ni(r,t)](r,t). 

Proceeding in a similar fashion, one can set up the Dyson-type integral equation for the 
third-order density response Eq. (lC9|) : 



n3 ^ = 6 J dy J dy> J ^ / '^ 3 H x >y>y^y / 0^s[ n l](y)^sK](y0^s[^l](/^ 

+ \l dy I dy ' I dZ I dz '^ { s ) ^ x ^y^y') U Ks[ni\{y)g X c{y\z,z')n l (z)n 1 (z') 
+ dy dy' dy"xf ] ' {x 1 y 1 y')U K s[ni]{y){w{y l ,y") + f xc (y' ,y"))n 2 (y") 
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+ \ I dy I dz [ dz' [ dz"x ( P{x,y)h xc (y, z, z' , z'^n^n^z^n^z") 



+ JdyJ dy' J dy" X ^(x,y)g xc (y,y , ,y")n 1 (y')n 2 (y") 
+ JdyJ dy'x^\x,y){w{z,z') + f xc {y,y'))n^y'), 



(Dll) 



where h xc is the third-order functional derivative of the time-dependent exchange-correlation 
potential with respect to the time-dependent densities 



h xc (y,z,z',z") 



S 3 U xc (y) 



5n(z)5n(z')Sn(z / ') 



(D12) 



"0 
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